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ABSTRACT 

We present a kinetic theory for the evolution of the phase-space distribution of dark 
matter particles in galaxy halos in the presence of a cosmological spectrum of fluctu- 
ations. This theory introduces a new way to model the formation and evolution of 
halos, which traditionally have been investigated by analytic gravitational infall mod- 
els or numerical N-body methods. Unlike the collisionless Boltzmann equation, our 
kinetic equation contains nonzero terms on the right-hand side arising from stochastic 
fluctuations in the gravitational potential due to substructures in the dark matter mass 
distribution. Using statistics for constrained Gaussian random fields in standard cos- 
mological models, we show that our kinetic equation to second-order in perturbation 
theory is of the Fokker-Planck form, with one scattering term representing drift and the 
other representing diffusion in velocity-space. The drift is radial, and the drift and dif- 
fusion coefficients depend only on positions and not velocities; our relaxation process in 
the quasilinear regime is therefore different from the standard two-body relaxation. We 
provide explicit expressions relating these coefficients to the linear power spectrum of 
mass fluctuations and present results for the currently favored cold dark matter model 
with a nonzero cosmological constant. Solutions to this kinetic equation will provide a 
complete description of the cold dark matter spatial and velocity distributions for the 
average halo during the early phases of galaxy halo formation. 



Subject headings: cosmology: theory — large-scale structure of universe — gravitation 
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1. Introduction 

Most mass in the universe is in the form of dark matter, and most dark matter is gravita- 
tionally clustered in the form of halos. Dark matter halos are therefore the building blocks of the 
universe. Knowledge of the formation and evolution history of dark matter halos is essential for 

an understanding of the halos' structure and dynamics, as well as the environment that harbors 
galaxies and clusters. Designs for dark matter search experiments also rely on the predicted spatial 
and velocity distributions of dark matter particles in our own galaxy. 

According to the theory of structure formation in current cosmological models, dark matter 
halos arise from tiny primordial fluctuations seeded in the matter and radiation fields in the early 
universe at, for example, the end of an inflationary era. These small inhomogeneities grow via gravi- 
tational instability, with the slightly overdense regions becoming denser and the slightly underdense 
regions becoming emptier. At a given time, there exists a spectrum of density fluctuations over 
a wide range of length and mass scales, whose spatial distribution is characterized by the power 
density in Fourier space P{k,t), the power spectrum. Dark matter halos therefore rarely form 
and evolve in isolation, growing instead by frequent accretion of smaller-mass halos and occasional 
major mergers with other halos of comparable mass. 

The earliest theoretical insight into cosmological halo formation came from the spherical radial 
infall model, which considers an isolated spherical overdense volume and describes the secondary 
infall of bound but initially expanding coUisionless matter onto this region (Gunn & Gott 1972). 
If the initial density perturbation spectrum is scale-free and the velocity distribution follows the 
Hubble law, the infall solution is found to approach a self-similar form in an Einstein-de Sitter 
universe. The resulting density profile asymptotically approaches p oc r~" with a > 2 (Fillmore 
& Goldreich 1984; Bcrtschinger 1985). The exact value of a depends on the index of the initial 
perturbation spectrum P{k) oc fc*^: it is isothermal (a = 2) for n < 1 and for n > 1 steepens to 
a K 3(n + 3)/(n + 4) (Hoffman & Shaham 1985) or 3(n + 3)/(n + 5) (Syer & White 1998; Nusser & 
Sheth 1999). The collapse of halos therefore retains memory of the initial conditions in the radial 
infall model. 

By contrast, modern cosmological N-body simulations find dark matter halos over a wide range 
of mass scales to follow a density profile that is nearly independent of the initial power spectrum 

and is shallower near the halo centers (with a ~ 1) than that predicted by the radial infall model 
(Navarro et al. 1996, 1997; Moore et al. 1998). Furthermore, the dark matter in the simulated halos 
is not entirely smoothly distributed spatially; instead, at least 10% of the halo mass is in the form of 
hundreds to thousands of smaller, dense satellite halos of varying mass (Ghigna et al. 1998; Klypin 
et al. 1999; Moore et al. 1999). Thus even though the analytic similarity solutions derived from one- 
dimensional infall onto a smooth overdensity have provided valuable insight into the evolution of 
self-gravitating systems, realistic halo formation is clearly more complicated. Numerical simulation 
has been the tool most frequently used in cosmology to study structure formation. It has provided 
powerful insight and occasionally surprising results (e.g. universal density profile) unexplainable by 
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the simulations themselves. Our understanding of structure formation can certainly be enhanced 

if a complementary tool is developed. Complementary approaches are also especially needed at 
this time to help us better understand and perhaps resolve the "cuspy halo" and "dwarf satellite" 
problems of the cold dark matter model. 

In this paper we offer a new perspective on the subject of dark matter halo formation by 
proposing an approach different from both the radial infall model and N-body simulations. We 
focus on the phase-space distribution of dark matter in galaxy halos instead of on the motion of a 
mass shell (as in infall models) or the orbit of a simulation particle (as in N-body). We attempt 
to understand the evolution of the phase-space distribution with a statistical description based on 
kinetic theory. Our method is to be contrasted with the standard method where the individual 
dark matter particles are described by the collisionless Boltzmann equation; N-body simulations 
are in essence a Monte Carlo method for solving this equation. Instead, we consider the phase 
space of dark matter particles in an average halo and develop a cosmological kinetic theory for halo 
evolution in the presence of a spectrum of density fluctuations. 

The essential features of our approach are represented by the fluctuation-dissipation theorem 
for time-dependent stochastic processes. To elaborate on this concept, we recall that stochastic 
modelling was first developed to explain the Brownian motion of pollen particles suspended in a 
fluid. (See § 5.3 for a more detailed discussion.) This phenomenon is now understood to arise from 
a fluctuating force on a Brownian particle due to many rapid collisions with the molecules in the 
surrounding fluid. These fluctuations cause successive small changes in the particle velocity, hence 
to a random walk. The consequence of many such random walks is diffusion, a form of dissipation. 

Due to the large number of molecules, it is impossible to solve the coupled equations of motion 
for the trajectories of all molecules and the pollen particle in the fluid. Fortunately, a stochastic 
treatment suffices. One mathematical prescription is provided by the Langevin equation (Langevin 
1908), which models the forces arising from the fluctuations as a stochastic term that is added to 
all other deterministic forces in the equation of motion for a particle. An alternative prescription 
is given by the Fokker-Planck equation (Fokker 1914; Planck 1917), which is an equation of motion 
for the probability density or distribution function of particles undergoing random walks. In this 
case, the collective changes in the velocities due to the fluctuations in the system are represented 
by a drift coefficient and a diffusion coefficient arising from the collision terms in the coUisional 
Boltzmann equation. 

Stochastic phenomena occur beyond molecular scales. In astrophysics, the best studied exam- 
ple is the evolution of dense stellar systems such as globular clusters (Chandrasekhar 1943; Spitzer 
1987; Binney & Tremaine 1988). A star in a globular cluster experiences a fluctuating force from 
many gravitational two-body encounters and suffers a succession of small velocity changes. The 
fluctuations in these systems arise from granularity in the mass distribution at the locations of 
the stars. The Fokker-Planck approximation is valid for these weak encounters. We explore this 
example in greater detail and contrast it with our model of galaxy halo evolution in § 5.4. 
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The key feature common to all of these systems is the presence of stochastic fluctuations that 
result in dissipation. In the cosmological setting addressed by this paper, the fluctuations arise from 
cosmological perturbations produced in the early universe. At early times, these perturbations are 
Gaussian and fully characterized by their power spectrum; at late times, the substructures in halos 
seen in recent simulations are the nonlinear counterparts of this spectrum of early ripples. These 
fluctuations scatter particles. If one considers an ensemble of halos, the masses and locations of the 
substructures (or their Gaussian predecessors) arc random variables. The substructures thereby 
introduce "noise" into the evolution of dark matter halos. As we show in this paper, the result is 
dissipation in the form of drift and velocity-space diffusion. 

Recent work has explored the use of the Fokker-Planck equation in cosmological problems. 
Evans & Collett (1997), for example, searched for steady-state solutions of the Fokker-Planck 
equation assuming the fluctuations arise from two-body encounters. It also assumed a simple 
isotropic phase space distribution dependent only on the total energy and a radial power law 
form for the potential and density. Under these assumptions, Evans k. Collett (1997) found that 
p oc r~^f^ was a steady-state solution under the effects of encounters. There was no dynamics in 
this work, however: the system was assumed to be in coUisionless equilibrium and the nature of the 
clumps producing the collisions was left arbitrary. It also remains to be determined whether galaxy 
halos, either dark matter or stellar, actually evolve to a = | cusp. Weinberg (2001) examined 
the dynamics of halos due to stochastic fluctuations from three different noise processes for the 
Fokker-Planck scattering term: the merging clumps were assumed to fly by with constant velocity, 
to follow decaying orbits due to dynamical friction, or to orbit inside halos. The results of this 
study were promising and supported the notion that the density profiles of halos were affected 
by mergers and substructures. The fluctuations in this work, however, were constructed by hand 
instead of being self-consistently calculated, and the results about a central cusp were inconclusive. 

Our objective here is to construct a kinetic equation in which the drift and diffusion terms are 
derived from realistic density fluctuations in current cosmological models. We will also retain the 
full phase space information by, in effect, using f{E,J,Vr,t) without averaging over the angular 
momentum or radial velocity. We also do not assume coUisionless equilibrium; instead we base our 
approach on exact dynamics. 

The organization of this paper is as follows: 

In § 2, we derive the kinetic equation relevant for the evolution of an average halo, starting with 
the coUisionless Boltzmann equation for individual dark matter particles. This procedure gives rise 
to a collision term on the right-hand side of the kinetic equation. This term represents a correlated 
force density and arises from fluctuations in the ensemble averaged gravitational potential due to 
clustering of the matter distribution. 

In § 3, we introduce the concepts and physical variables needed to apply this equation to cos- 
mological problems. The key step is in relating the phase space density f{f,v,t) to the probability 
densities of two standard variables in cosmology: the mass density p and velocity v. We then focus 
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on the quasi-linear regime in which the fluctuations about the mean density and Hubble flow of 

an isotropic and homogeneous universe are small. This focus allows us to relate p and v to the 
density and displacement flelds 8 and V', which are Gaussian random fields in standard cosmological 
models. 

We leave the most mathematical part of the paper to Appendix A, where we describe how 
to use the statistics of constrained Gaussian random fields to express the probability densities of 
§ 3 in terms of the (constrained) means and covariances of the density and displacement fields. 
The idea of a constrained field is central to this study because our aim is to construct a kinetic 
theory for dark matter particles that will later reside in a collapsed halo rather than at any random 
location in the universe. The relevant perturbations are for a density field constrained to have 
certain properties, e.g., a specified height or gradient at some point in space. The spectrum for 
such constrained fields is related but not identical to the familiar linear power spectrum P{k) for 
unconstrained Gaussian density fluctuations. 

In § 4, wc put the final expression for the correlated force derived in Appendix A back into 
the kinetic equation, and show that the resulting equation is of the Fokker-Planck form where one 
term represents a drift force and the other term represents diffusion in velocity-space. We work out 
these diffusion coefficients for the currently-favored cold dark matter plus a cosmological constant 
(ACDM) model, and provide explicit expressions for these coefficients in terms of the linear matter 
fluctuation power spectrum. We also present the results for the coefficients computed for this 
model. 

§ 5 discusses further the physical meaning of the results in § 4. Wc interpret the diffusion terms 
by taking velocity moments of the kinetic equation. To place our cosmological kinetic equation in 
a broader perspective, we discuss in more detail the classical examples of Brownian motion and 
globular clusters, and point out the similarities and differences among the three cases. § 6 provides 
a summary and conclusions. In Appendix B we investigate the relaxation to stationary solutions 
of our cosmological Fokker-Planck equation. 



2. Derivation of the Kinetic Equation 

We derive the kinetic equation for galaxy halo evolution following the methods presented 
by Bertschingcr (1993) for the derivation of the BBGKY hierarchy. The starting point is the 
one-particle phase space density for dark matter particles in a single halo, the Klimontovich (or 
Klimontovich-Dupree) density (Klimontovich 1967; Dupree 1967) 

fK{r,v,t) = m'^SBir - fi{t)]dB[v - Vi{t)] (1) 

i 

where (5d is the Dirac delta function and the subscript K denotes the use of the Klimontovich 
density. We suppose that each halo is made of equal-mass particles of mass m, for example, the 
elementary particles of non-baryonic dark matter models. For convenience, we use velocity rather 
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than momentum and we choose units so that fd^v is a mass density. We use proper coordinates 
except for § 3.2 and Appendix A, where we will switch to the appropriate cosmological (comoving) 
coordinates. 

The Klimontovich density is simply a way of casting the exact trajectories of all N particles 
in a system into the form of a phase space density. The number of particles in d^rd^v about 
(r, v ) \s fKd^r d^v, which equals either or 1. The Klimontovich density consists of delta functions 
because, when one is dealing with a single system (as opposed to an ensemble), every particle has 
a well-defined position and velocity. The reason for introducing this quantity is that it gives a 
particularly clean way of going from exact N-body dynamics to a statistical description by taking 
expectation values of /k or products of more than one /k over an ensemble of N-body systems. 
This approach to kinetic theory is different from the usual method that begins with the N-particle 
distribution for an ensemble. 

The Klimontovich density obeys the Klimontovich-Dupree equation 

dfK ^ ^ 5/k 9/k „ 

^ + -^+5K-^=0, (2) 

where 

Mr,t) = -Gm^^^ = -G J d'w' Mw' ,t)^^^ . (3) 

For notational convenience we have grouped together all 6 phase space coordinates into w = {r, v}. 
Eq. (2) is exact, as may be verified by substituting Eq. (1) and using the equations of motion for 
individual particles, 

df ^ dv _ 

Tt=''^ -dt=^- 

In practice, the gravitational force may be softened by modifying Eq. (3) to reduce artificial two- 
body relaxation, as is done in N-body simulations. Eqs. (1) and (2) are entirely equivalent to Eqs. (4) 
for one system. The Klimontovich phase space density retains all information about a particular 
halo because it specifies the trajectories of all particles. Eqs. (l)-(4) are formally equivalent to a 
perfect N-body simulation. That is not what we want. Rather than giving a perfect description of 
a single halo, we average over halos to obtain a statistical description of halo evolution. While it is 
impossible to completely describe a single galactic halo of dark matter particles, by using statistical 
mechanics methods we can describe the average of infinitely many halos! 

Before proceeding further, we make one important modification to the usual procedure in 
kinetic theory. Eqs. (4) assume that we work in an inertial frame. This is satisfactory for an 
isolated halo but not for a halo formed in a cosmological context where neighbors cause it to orbit 
in the net gravitational field produced by all matter. The center of mass of a halo typically moves 
much farther in a Hubble time than the size of the halo itself, even in an inertial frame in which the 
halo begins at rest. For example, for the cold dark matter family of models, the displacement of a 
mass element from its initial position (in comoving coordinates) is typically several Mpc today. If 
we statistically average over an ensemble of halos each of which wanders this far in a Hubble time. 
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the phase space density will be artificially smeared out, making the resulting average "halo" much 
larger than any realistic halo. We avoid this by doing the same thing that simulators do when they 
construct halo profiles: we work in the center-of-mass frame of each halo. 

Indeed, at each moment the halo center-of-mass must be at rest in our reference frame. This 
frame is non-inertial because of the gravitational forces acting on the halo. Thus we must modify 
the kinetic equation to apply to an accelerating frame. Fortunately, the accelerations are small and 
the frame is extended only a few Mpc, so we do not have to worry about relativistic effects. We 
simply include the inertial ( "fictitious" ) forces associated with working in an accelerating frame. 

We choose to define r = as the center-of-mass of each halo. Transforming Eqs. (4) to the 
center-of-mass frame is very simple: we replace g by the tidal gravitational field relative to the 
center of each halo, g{f) — g{0), and then subtract the center-of-mass velocity from the initial 
velocity of each particle. Eq. (2) is unchanged except that g is replaced by 

gMr,t) = -GmY: + ^) = -G j d'w'Mw',t) + ^) . (5) 

The subscript T is a reminder that this is the tidal gravitational field. 

Now we return to the usual procedure for deriving a kinetic equation. The Klimontovich 
density is the phase space density for a single realization of a halo. We might imagine performing 
many N-body simulations of different halos, each with slightly different initial conditions. If we 
average the Klimontovich density over this ensemble of halos, the result (at least, in the limit where 
the size of the ensemble becomes infinite) is no longer a sum of delta functions but instead is a 
smooth one-particle distribution function: 

f{w,t) = {Mw,t)), (6) 

where angle brackets denote the ensemble average. Recall that for notational convenience we have 
grouped all six phase space variables into w. 

Our goal is to derive a tractable kinetic equation for /. The obvious next step is to take 
the ensemble average of Eq. (2). To proceed, however, we will need the two-particle distribution 
function f2{wi,W2,t), which follows from averaging the product of Klimontovich densities at two 
points: 

(/K(w^i,i)/K(w^2,i)) = <5d(i(^i - iy2)/(tyi,i) + /2(w^i,'i«2,i) • (7) 

The delta-function term 5d{wi — W2) = SD{ri — r2)ST>{vi — V2) corresponds to the case in which the 
two points lie on top of one particle. Two distinct particles give the /2 term, which may always be 
written as a product term plus a two-particle correlation: 

f2{wi,W2,t) = f{wi,t)f{w2,t) + f2c{wi,W2,t) . (8) 

This equation defines the two-point correlation function in phase space, f2c- 
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The averaging we perform is similar to the averaging used in defining ordinary clustering 

correlation functions in cosmology. The differences here are that we are not averaging over arbitrary 
systems but only over those that lead to formation of a halo at r = 0, and that we retain velocity 
information. Had we averaged over all space and integrated over velocities, we would have gotten 
/ d^v f = p and / Svi cfv2 /2c = P^S, where p is the cosmic mean mass density and ^ is the usual 
spatial two-point correlation function. 

Now, by taking the ensemble average of Eq. (2) we arrive at a kinetic equation for the evolution 
of the average halo phase space density (Bertschinger 1996): 

— + V ■ ^ + qt ■ -rz = ^ - F , (9) 

dt df dv dv ' 

where we define the ensemble-averaged tidal field 

gTivM) ^ iffKTim =-gJ ^w'f{w',t) + ^) (10) 

and the correlated force density 

F{w,t) = Cov[gKTir,t),fKiw,t)] = -G j (fw'f2c{w,w',t) (^ |!j!,|3 + ^) • (H) 

The covariance of two random variables is denoted Cov[A, B] = {{A — {A)){B — {B))) = {AB) — 
{A){B). The product of ^kt and /k involves a product of two Klimontovich densities, for which 
we have used Eqs. (7) and (8). The one-particle discreteness term in Eq. (7) makes no contribution 
because the particle has no self-force. 

Eq. (9) is the first BBGKY hierarchy equation (Ichimaru 1973; Davis &: Peebles 1977; Peebles 
1980). In cosmological applications this equation is usually derived starting from the Liouville 
equation for the N-particle distribution function. In order to clarify the averaging that is done 
here, we have started instead from the one-particle Klimontovich density. The result, Eq. (9), is an 
evolution equation that looks very much like Eq. (2). By averaging over halos, however, we have 
introduced a correlation integral term on the right-hand side. The average halo therefore does not 
evolve according to the Vlasov equation; instead, it evolves according to a kinetic equation with an 
effective collision term. Eqs. (9)-(ll) are exact. No approximation has been made yet. 

Eq. (9) is a continuity equation for density in phase space. Without the right-hand side, it 
says that phase space density is conserved along trajectories given by Eq. (4) modified for the tidal 
force. The right-hand side gives a correction due to the fact that fluctuations about the ensemble 
average lead to correlated forces. As we will see, the effects of these fluctuations lead to dissipation. 
Because the correction term is a divergence, the mass density p{r, t) = J f (Pv is not dissipated but 
instead is conserved. 



In most applications of kinetic theory, the right-hand side of Eq. (9) arises from particle 
collisions. In the usual treatment of a gas of particles with short-range forces, one assumes /2c = 
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always except during particle collisions. This is the assumption of molecular chaos that Boltzmann 
introduced to derive his celebrated kinetic equation. Our situation, however, is different: we cannot 
neglect /2c, because it can be much larger than the single-particle term in Eq. (8) due to strong 
gravitational clustering. 

Hcuristically, /2c describes the substructure within a galaxy halo at the two-point level. (Higher 
order correlation functions would be needed for a complete description.) Current cosmological mod- 
els and observations imply that the initial density field has fluctuations that cause the formation 
of many small halos which subsequently merge into present-day halos. The lumpiness of the mat- 
ter distribution represents a fluctuation about the ensemble average density field. Through the 
correlated force density F, these fluctuations cause changes in the energy and angular momentum 
of individual particle orbits that are important for the evolution of the one-particle distribution 
function. 

Eq. (9) is useful only if wc find an expression for /2c. One way to proceed would be to derive a 
kinetic equation for /2c by taking the ensemble average of Eq. (2) multiplied by the Klimontovich 
density. This leads to the second BBGKY equation, which depends on the three-point correlation 
in phase space, which obeys the third BBGKY equation, and so on. Thus we either have an infinite 
chain of coupled equations of increasing dimensionality, or we must close the system. Davis & 
Peebles (1977) handled this problem in their study of gravitational clustering by writing the three- 
point function in terms of products of two-point functions. We will proceed differently, by seeking 
to close the hierarchy at the first level, Eq. (9). 

To summarize this section, wc have derived a formally exact evolution equation for the average 
halo. The one-particle phase space density is conserved up to a fluctuating force term arising from 
the two-point correlation function in phase space. Physically this term represents the effects of 
substructure within and around halos, which scatter particles to new orbits. In the next two 
sections, we will show that we can in fact express /2c in terms of the one-particle phase space 
density / using a relation that is exact to second order in cosmological perturbation theory. This 
relation will enable us to close the first BBGKY hierarchy. 

3. Distribution Functions from Probability Theory 

Solving our evolution equation in Eqs. (9)-(ll) requires specifying the phase space density / 
at an initial time and the two-particle correlation /2c at all times. In this section we show how 
to compute these quantities from the probability distributions of mass density and velocity at one 
and two points in space. 
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3.1. Probability Density vs. Phase-Space Distribution Function 

We wish to derive expressions that would relate f(wi,t) and f2c(,wi,W2,t) to probability dis- 
tributions of mass density and velocity, because these latter quantities are naturally specified by a 
cosmological model for structure formation. The method is to examine the mass and momentum 
in a small vohimc of space and to construct ensemble averages of the Klimontovich density using 
the probability distribution of density and velocity. 

We use the symbol p{x, y,---) to refer to the probability density with respect to its arguments, 
where the joint probability distribution for ...) is dP{x,y, ...) = p{x,y, . . .) dx dy ■ ■ ■ . We 

assume that these probability distributions are always normalized to unit total probability. The 
variables relevant for this paper are the densities p and velocities v at two points fi,f2 and at 
one time t: pi = p{fi,t), vi = v{fi,t), p2 ^ p{f2,t), and V2 = v{f2,t). Initially we assume that 
the velocity field in a given realization is single-valued in space. We use the following probability 
distributions: 

dP{pi) = pipi)dpi 
dP{pi,vi) = p{pi,vi)dpid.^vi 



dP{pi,vi,p2) = p{pi,vi,p2)dpid'^vidp: 



'2 

dP{pi,Vi,p2,V2) = p{pi,Vi,p2,V2)dpid^Vidp2d^V2 , (12) 

where the different probability densities are related to one another by 

PiPl,Vl,P2) = J d^V2pipi,Vi,p2,V2) 
P{Pl,Vl) = J dp2Pipi,Vl,P2) 

P{pi) = J d^vip{pi,vi) 

1 = J dpipipi). (13) 

These probability distributions depend on time but we suppress the f-dependence for clarity. 

For a small volume in phase space, f{wi,t) d^r\d^vi equals the mean mass contained in 
d^ri d^vi. The same quantity can be calculated by averaging pi d^ri over density at fixed velocity 
using dP{pi,vi), i.e., J pid^r\dP{p\,vi). Equating the results gives 

/*00 /"OO 

f{wi,t)= dpipip{pi,vi) = dpipip{pi\vi)p{vi) = {pi\vi) p{vi) . (14) 
Jo Jo 

We have used the conditional probability density p{pi\vi) to define the conditional mean {pi\vi). 
Note that the probability densities depend on (ri,t) through pi = p(fi,t) and vi = v{fi,t), so 
Eq. (14) yields the desired one-particle phase space dependence on wi = {fi,vi}. 



This argument is easily extended to the two-point phase space density, using the joint proba- 
bihty distribution of density and velocity at two points. As Eq. (11) shows, we need f2ciwi,'r2,t), 
which depends on vi but not V2, i.e., we need the velocity at only one point. Recalling that 
f2ci'Wi,W2, t) = f2i'Wi,W2,t) - f(wi,t)f{w2,t), we obtain 



Although we have derived Eqs. (14) and (15) assuming a single-valued velocity field in each re- 
alization, they are also valid when there is a distribution of velocities at each point. One simply 
interprets pi as the total mass density while vi is a single particle velocity. Therefore these equations 
are fully general. 



At high redshift, before dark matter halos formed, the density and velocity fields were only 
slightly perturbed from a uniformly expanding cosmological model. It is thought that the fluc- 
tuations of density and velocity responsible for all cosmic structure were Gaussian random fields 
produced during the inflationary era of the very early universe. These fluctuations provide us with 
both the initial conditions and the random element leading to a probabilistic description for halo 
formation. Our goal now is to obtain expressions for Eqs. (14) and (15) at early times while the 
density fluctuations are small. 

Before proceeding further, we review the description of density and velocity fields at high red- 
shift in an almost uniformly expanding cosmological model. First, the mean expansion is described 
by the cosmic expansion scale factor a{t) which is related to the redshift by 1 + z = l/a{t). The 
Hubble expansion rate is H(t) = d/a. The matter contributes a fraction i^mit) to the closure 
density. We account for the mean expansion by defining new coordinates conventionally called 
comoving coordinates: x = r/a{t). 

Adopting the standard cosmological model, we assume that the dark matter is cold. This 
assumption means that the velocity field is single-valued in space at early times before halos form. 
We can then write the density field p and velocity field v in terms of perturbations 5 and ip around 
a uniformly expanding background: 



Here h{t) = a{dlnD/d\na) ^ a[nm{t)]'^'^ is the growth rate for density fluctuations with time- 




(15) 



3.2. Cosmological Variables 



p{r,t) = p{t) [l + 6{x,t)] 



v{f, t) = H{t) r + b{t)^f{x, t) 



(16) 
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dependence S{x,t) oc D{a). We emphasize that Eq. (16) does not assume that the density fluctua- 
tions are small; it only assumes that the velocity field is single- valued. We are implicitly assuming 
that all of the matter may be described as a single cold non-interacting fluid (i.e. we assume that 
baryons move like cold dark matter over the length scales of interest). 

The computation of the functions a{t), H{t), p{t), D{t), and b(t) is standard in cosmology and 
will not be described here. At early times when 5^ <C 1, mass conservation relates the density and 
velocity perturbation flelds by S = —{d/dx) ■ ip, or 

^>"'" = -/^|^*(^-'"- (1^) 

Now we make a crucial step by changing variables from (f,v,p) used previously to (x,^,6). 
The reason for doing so is that the description in comoving coordinates x and perturbation variables 
tp and S is much simpler in the small perturbation limit. Following the notation in § 3.1, we write 
5i = 5{xi,t), 82 = 6{x2,t), and ^pl = tlj{xi,t). These variables have normalized joint probability 
densities p{Si), p{Si,tjji), and pi^i, 1^1,62) in analogy with Eqs. (12) and (13). 

To change variables from {v, p) to (V', 5), we use dp = pdS, d^v = (Hb)^ d^ip, and dp d^vp{p, v) = 
d6 d^ijjp{6,il^ ). Eq. (14) then becomes 

fiwi,t) =p[l + {Si\ifi)] p{ifi) (Hb)-^ , (18) 

where 



{6i\^i) = J d6ip{6i\iJi)5i , p{Si\i;i)=p{5i,^i)/p{^i) , (19) 
and p{Si\'tpi) is the conditional distribution of Si. Similarly, Eq. (15) becomes 

f2c{wuX2,t) = p2 \{S2\i^,) - {62) + ((^K^slV^l) - {SMl){S2)] p{i^l){Hb)-^ . (20) 



Although the derivation of Eqs. (18) and (20) assumed that the velocity fleld is single- valued, 
they are valid in general provided that one interprets Si and S2 as total mass density perturbations 
(not necessarily small) while tpi is proportional to the single particle velocity. Eqs. (18) and (20) are 
equivalent to Eqs. (14) and (15) expressed in different variables. Their utility becomes evident in 
Appendix A where we evaluate the probability distributions of density and velocity in cosmological 
models. 

All that remains is to calculate the expectation values in Eqs. (18) and (20) and to substitute 
the latter into Eq. (11). This calculation is somewhat complicated and is presented in Appendix A. 
Part of the complication arises because of an important detail we have so far neglected: halos form 
preferentially in regions that are initially overdense. Our calculations assume that a halo forms at 
r = 0. Consequently we should not evaluate the expectation values for arbitrary points in space 
but rather for those that satisfy appropriate constraints. 
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Strictly speaking it is impossible to specify in advance the regions that will form halos in a 
cosmological model with random initial fluctuations because the dynamics of gravitational cluster- 
ing is too complicated. A simple model, however, has been found to give reasonable agreement 
with numerical simulations: halos form at maxima of the smoothed initial density field (Bardeen 
et al. 1986). These constraints are found to modify the expectation values in Eqs. (18) and (20). 
To emphasize the application of constraints, we rewrite these expectation values with a subscript 
C: 



and 



f{wi,t)=p l + (5i|V'i)c v{^i){Hh) 



f2c{wi,X2,t) = (^2|V'i)c - {62)0 + {h52\tp\)c - {h\'^i)c{h)c Pitpi){Hb) 



(21) 



(22) 



Eqs. (21) and (22) are the main results of this section. We will use them to specify the initial 
conditions on / and to determine the correlated force density in Eq. (11). 

The next step is to evaluate the necessary constrained expectation values using the standard 

cosmological model. We present the details of this calculation in Appendix A, where we choose to 
err on the side of simplicity by following Bardeen et al. (1986) and using only the zeroth and first 
derivatives of the smoothed initial density field A(x ) as constraints. Detailed comparison has shown 
that this simple model does not predict perfectly all halo formation sites (Katz, Quinn, & Gelb 
1993), and more complicated constraints work better, e.g. Monaco et al. (2002). The procedure of 
Appendix A can presumably be extended to incorporate more complicated constraints. 



4. Fokker-Planck Equation 

After having devoted § 3 and Appendix A to probability and statistics we are now ready to 
return to the kinetic equation in § 2 and write down an explicit expression for the correlated force 
density F{w, t) on the right-hand-side of the kinetic equation (9). The approach we have followed is 
exact rather than phenomenological — our starting point was the first BBGKY hierarchy equation, 
an exact kinetic equation. The only approximation we have made is to assume that the density 
and velocity fields are Gaussian random fields. The resulting two-particle correlation function, 
Eq. (A22), depends in a simple way on the one-particle distribution p(V'i)- As we will show next, 
this simple dependence leads to the Fokker-Planck equation. This is a remarkable and non-trivial 
result. The Fokker-Planck equation is a great simplification of the exact — and intractable — 
BBGKY hierarchy. Usually the Fokker-Planck equation is assumed to hold a priori, and then its 
diffusion coefficients are estimated. Here, on the other hand, we have derived the Fokker-Planck 
equation. Along with it we obtain the diffusion coefficients in the quasi-linear regime, as follows. 
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4.1. Diffusion Coefficients 

Substituting Eq. (A22) into Eq. (11) using Eq. (15), converting to comoving coordinates, using 
Eq. (17) to relate to v, and using Eq. (A21), to third order in S we get 

F{w, t) = -ATTGpa { [C{6, - C{5, ■ C-\'^,^) ■ C{^, ^o)] f 

+ (Hb) [C(V^, ) - Ci^Po, i> ) - (V^)c ® C{5, v;)] ■ + 0(6^) , (23) 

where the covariancc matrices C of variables 6 and ip are given in Eqs. (A14) and (A15). Here, 
i/j = ip{x,t) and tfjQ = ■0(0, t). Note that all the terms with subscript arise from our use of the 
tidal field rather than the total gravitational acceleration. 

Substituting Eq. (23) into Eq. (9), in proper coordinates our kinetic equation in the quasilinear 
regime reduces to the Fokker-Planck equation 

dt ^ dr dv dv \ dv J 

specified by a vector field A and a tensor field D. Together, A and D are called diffusion coefficients. 
Heuristically, A represents a drift force {A has units of acceleration) while D represents the diffusive 

effects of fluctuating forces. The flux density in velocity space is Af — D • df /dv. Keeping only 
terms up to second order in the density perturbations, these quantities are given by 

A{r, t) ^ -AirGpa [c{5, ^o) - C{d, ) • C-^tf, ^J) ■ C(0, ) 

= -CoYi5,go\^) = CoYi5,gT\v) (25) 

and 

B{r,t) ^ 47:GpaHbC{ip - i^o,4^) = Cov{gT,v) . (26) 

In Eqs. (25) and (26) we have also used the fact that in the quasilinear regime the peculiar gravity 
and displacement are related by ^ = AirGpa^ (Zel'dovich 1970), and the tidal field is = g{f, t) — 
go, where = 9(0, t). In writing the second line of Eq. (25) we have used Eq. (A7) with Yb = {5, tpo} 
and Ya = ■0— {t/j)c- Thus, Cov{6,go\v) should be understood to be subject to the same extremum 
or other constraints as e.g. the constrained covariancc function C{S,tp). The displacement and 
velocity are related by the second of Eqs. (16). Thus, in linear theory g v so that conditioning 
on V makes the fluctuation in g vanish, allowing us to replace —Cov{6, go\^) by Cov(5, ^xl^'") in 
Eq. (25). In all cases, the covariances are taken subject to the extremum constraints discussed in 
Appendix A. 

It is worth noting that had we used the total gravitational acceleration g{f, t) instead of the 
tidal field gT{f,t), the drift vector A would have been zero and D = Cov{g,v) would have been 
larger. We conclude that tidal effects are crucial for a proper description of halo evolution (Dekel, 
Devor, Sz Hetzroni 2003). This is evident in the final expressions for A and D. After the lengthy 
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derivations, these expressions are remarkably simple, suggesting that they may have greater validity 
than the case of small-amplitude Gaussian fluctuations that we have analyzed here. 

In going from Eqs. (A21) and (23) to Eqs. (25) and (26) we have dropped the contributions 

made by {5)c and {ip )c because they do not contribute at second order in 5. At this order, A and 
D are independent of the value of the smoothed density Aq used for the constraint. However, the 
constraint enters in third order. 

Carrying out the derivation to one higher order in 5 requires using the following result for 
Gaussian random variables: Cov(^i?,C) = {A)Cov{B ,C) + {B)Coy{A,C). Applying this with 
S = (1 + 5)~^ = 1 — 6 + 0{6'^), after some algebra we find 



l + <5 



A{r, t) = Gov <5, V + 0(5^) , D(r, t) = Gov rf^,,v + 0((5^) . (27) 
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1 + 5' 



Although these results are correct for Gaussian fluctuations to third order in the density pertur- 
bations, nonlinear evolution will generate non-vanishing irreducible three-point correlations (i.e., 
non-Gaussianity) that have been neglected here. For example, nonlinear evolution generates cor- 
rections to the Zel'dovich approximation so that velocity and gravity are no longer proportional to 
each other as we have assumed. For this reason, Eqs. (27) should not be regarded as being more 
accurate than Eqs. (25) and (26). 

Despite their limitations, Eqs. (27) are useful in showing the qualitative effect of applying an 
initial constraint. Choosing the proto-halo to have Aq > means that typically > in the 
proto-halo, which decreases A and D compared with Eqs. (25) and (26). Conversely, the diffusion 
coefficients are enhanced in void regions. Note also that the constraint value appears explicitly in 
the initial conditions for / in the linear regime through the term {5)c iu Eq. (A21). 

Eqs. (25) and (26) imply that, in the quasilincar regime, the drift and diffusivity are inde- 
pendent of the velocity within the halo. (Although A is equivalent to a covariance at fixed it is 
independent of the particular value of v.) This follows because /2c depends on V'l = {v\ —Hxi ) / (Hb) 
only through p{ijji) and dp/dijji in Eq. (A22). As noted above, at second order in perturbation 
theory the drift and diffusivity are also independent of the initial overdensity of the peak Aq (but 
they depend on the smoothing scale used to find extrema). Thus, both the Fokker-Planck equation 
and the diffusion coefficients appear to be remarkably robust quantities. 

To make further use of Eqs. (25) and (26), we need to write out the coefficients explicitly. 
Using Eqs. (A14), (25), and (26) we first factor out the dependence on cosmology by writing 

A{r,t) = 4TrGpaAr{r,t)f, (28a) 
D(r,t) = 4TrGpaHb[Dr{r,t)r^r + Dt{r,t){I-r^r)], (28b) 
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where Dr and Df are scaled radial and tangential diffusion coefficients. The scaled variables are 

Ar{r,t) = -rr]{r) - ^ 

^1 

d'y{r) 2 d'y{r) \rf]{r) 



,7(0) d{(r) ^ J^dr^ |«) ^ J^dgr) ^^^^^ _ ^^^^^ 



erf dr 



Cr{r) dr 



dr 



Dr{r,t) 
Dtir,t) 



Cr{r) 

ct{r) ■ 



dr 
7(r) 



a, 



d'j{r) 
7(r) 



77(7^ 



2 



■^(r) - 2;/(r) 



(29b) 
(29c) 



Here, r is the comoving distance from the center of the halo. The time dependence arises only 
because the functions appearing in Eqs. (29) depend on the time-dependent power spectrum of 
density fluctuations. These functions are summarized here for convenience: 



ry(r) = 
m = 



-^P(fc)^ 



(27r)3 ^ kr ' J (27r)3 ^ P ' 

^ Pik)WR{k)Mkr) , fi{r) = j -^P{k)WR{k)^-^ 
d^k „._„,o.,. o 1 r d^k 



kr 



(30a) 
(30b) 



where we recall that (Jq and af are the covariances of the constraints Aq and Aq defined in Eqs. (A8), 
and WR{k) is a smoothing window defined in Eqs. (A3) and (A4). We use the Gaussian window 
WR{k) = exp(-A;2i?2/2). 

To gain some insight into the significance of these results, let us consider a proto-halo at 
high redshift. If we suppose that the gravity field is dominated by its small-scale components, 
we may replace by g. If we further ignore the effects of initial constraints, in linear theory 
Dr Df ^ \Ha1 where is the one-dimensional velocity dispersion. In one realization, the forces 
are deterministic, but taken across the ensemble they are random and fluctuating. A diffusivity 
D ~ Ha^ implies that these fluctuations cause the typical particle to change its velocity by order 
itself in about a Hubble time. There is also a net drift force A with a dissipative timescale of 
about dynamical times, where (5 is a typical magnitude of the density fluctuations. As a result, 
friction and diffusion will significantly alter halos within a Hubble time of their collapse at high 
redshift. We conclude that it is not valid to model halos using the spherically symmetric Vlasov 
equation. This conclusion has also been reached previously by Antonuccio-Delogu k, Colafrancesco 
(1994). 

To summarize, Eqs. (24)-(26), and their more detailed version, Eqs. (28)-(30), are the main 
results of this subsection. Eq. (24) has the standard form of a Fokker-Planck equation. The Fokkcr- 
Planck form arises naturally and Eqs. (25)-(26) are exact for Gaussian fluctuations up to second 
order in perturbation theory. The diflFusion coeflficients A and D are due to fluctuations in the 
gravitational field over the ensemble of halos. The drift vector A represents a source of gravitational 
acceleration in addition to the acceleration produced by the mean matter distribution. The 
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diffusion term — D • df/dv causes the temperature (i.e. velocity dispersion) to increase (for a 
positive definite D). 

4.2. Results for ACDM Model 

Having derived analytic expressions for the diffusion coefficients in the previous subsection, 
we now evaluate them for the currently favored cosmological model at high redshift. Figures 1 
and 2 show A^, D^, and Dt as a function of radius r (measured from the halo center) for a range 
of Gaussian smoothing lengths R in Wr(/c). We obtain these results by numerically integrating 
Eqs. (29) and (30) with the linear power spectrum P{k) for a fiat ACDM cosmological model with 
{nm,nA,nb,h) = (0.3,0.7,0.05,0.7). 

The figures show a number of interesting features. First, the diffusion coefficients are negative 
over some regions of radius for some parameters. As we will show below, A is the gravitational 
acceleration on particles in the halo caused by drift. If the friction were produced by the Chan- 
drasekhar mechanism of tidal wakes, one would expect A \\ —v (see § 5.4). However, here A is 
independent of v and the particle mass. The source of friction is not the tidal wake of a flow past 
a moving massive particle. Instead it is the fluctuations of the gravitational field of a Gaussian 
random process, and these fluctuations can lead to a radially outward or inward force. We discuss 
this further in § 5. 

More surprising at first is the fact that the diffusivities Dr and Dt can be negative. Negative 
diffusivity causes the velocity dispersion (or temperature) to decrease, indicating an instability. 
Such behavior is possible for gravitational systems. In fact, this result is to be expected in the 
quasilinear regime of gravitational instability because second-order perturbations enhance gravi- 
tational collapse of dark matter halos (Peebles 1980; Jain & Bertschinger 1994). In the strongly 
nonlinear regime, after halo virialization, we expect the diffusivities to become positive. 

To elucidate the origin of the negative diffusivity, we show in Figures 3 and 4 that once the 
extremum constraint Aq = VA(0) = is removed (by setting aj"^ = 0), we obtain < and 
Dr,Dt > for all radii. The extremum constraint is therefore crucial to generating positive Ar 
and negative diffusivity seen in Figures 1 and 2. This effect is easy to understand: if one stands 
at a random place in a Gaussian field, nonlinear effects can speed up or slow down gravitational 
collapse. Gravitational instability, however, speeds up the evolution in the vicinity of a density 
extremum. 

Another obvious systematic effect with the diffusion coefficients in Figures 1 and 2 is their 
dependence on the smoothing radius R used to set the initial constraints. For the power spectrum 
assumed here, the maximum of A^ occurs at about 2R. In the limit -R ^ 0, the variances ctq and 
(jf diverge, causing the constraint terms in the diffusion coefficients (i.e., those terms proportional 
to and cr]~^) to vanish. 
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Figure 2 shows that the diffusivity is isotropic at the halo center but shows radial anisotropy 
outside the center. At small radii the tangential diffusivity is larger in absolute value than the 
radial diffusivity. Tangential diffusion gives the infalling matter angular momentum. Our treatment 
therefore provides a new statistical description of the growth of angular momentum by tidal torques 
(White 1984) and may enable in principle an explanation for the universal angular momentum 
profile found recently by Bullock et al. (2001). 



4.3. Asymptotic Behavior 



It is instructive to work out analytically the asymptotic behavior of the functions in Eqs. (29) 
for small and large r. For small r we find, 

"3^2(0) 



Ar{r,t) = 

Dr{r,t) = 

Dtir,t) = 
where we have defined 



afj{0) + 



1 



^1 

m 



-\ 2 



3^-1 + l^m 
4-.^ + l-m 



m 



+ 0(r^ 




+ 0(r^ 



i2n) 

m = 3f?(0) = 



3 m , 



a = 



J(fkk'^P{k)WR{k) 
J (Pkk'^P{k)Wl{k) 



(31a) 
(31b) 
(31c) 

(32a) 
(32b) 
(32c) 



Eqs. (31) show that if is finite, then as r ^ 0, drops to zero while Dj. and Dt approach a 
negative constant as discussed above, — f/^(0)/(Tf . This central value is decreased as the smoothing 
scale is decreased (i.e. higher do and ai) as shown in Figure 2. As i? — 0, we have gq as and 
o"i — GO as long as din P/ dink > — 5 as A; — oo. In the limit o"i — > oo (no extremum constraint), 
Dj. and Df are non- negative everywhere (cf. Figs. 2 and 4). Thus, negative diffusivity (instability) 
is produced by constraining the slope of the density field. Constraints have the opposite effect on 
Ar'. as the comparison of Figures 1 and 3 confirms, the contributions to Ar from constraint terms 
arc positive. 



For large r, both Cr{r) and Ct{r) approach cr^, and 7(r) oc 1/r (assuming P on k as k 



Dr{r,t) 



J2 



Dt{r,t) 



0), so 
(33) 



Figure 2 shows that Dr rises slightly above cr^ at a few hundred Mpc (for the standard model power 
spectrum) before reaching (tL This is because = Cr — d'y/dr and dj/dr is negative for this range 
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of r. The drift Arir^t) approaches zero at large r. In practice, the large-r hmit is unimportant 
because halo infall occurs only from within about 10 Mpc. 



4.4. Power-law Models 



Another way to gain insight into the diffusion coefficients is to examine their behavior for 
scale-free spectra P{k) oc k^. Many of the integrals are analytic in this case. 



We define a set of functions 



Jo 



ji{bx)exp 



1 



X ] dx . 



(34) 



The integral converges for a + I > 0. Using properties of the spherical Bessel functions ji{x), we 
obtain the following useful relations: 



bGi+i{a,b) = {a + l-l)Gi{a-l,b)-Gi{a + l,b) , 
bGi^i{a,b) = {2 + l-a)Gi{a-l,b) + Gi{a + l,b) , 
Go(a,0) = 2(«-2)/2r 



We also use the following relation, valid for a + / > 0: 
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(a)= lim [ x''~'^jl{x)exp(-^^dx = 2''-'^^/^T 



a + l 



3 + l-a 



(35a) 
(35b) 
(35c) 



(36) 



For Gaussian-tapered power-law spectrum P{k) = 27r^BA;" exp(— /c^i^g) and Gaussian smooth- 
ing Wnik) = exp(-A;2i?2/2), Eqs. (30) give, for i?o ^ and n > -3, 



rr]{r) 
7(r) 

d'f 
dr 

7 



Bffi(n + 2)r-("+2) , 
Bgi{n)r~'^ for n > — 1 , 

^n5i(n)r-("+i) + ir(^) R-^-+''> 



B 



J2 



= B 



n/-l : 
n = —1 

n / -1 
n = — 1 



4 



\ + \C + \ log(r/i?o) , 

5^(,,),-(n+l)+lr(H±l) 

I + iC + i log(r/i?o) , 
BGo{n + 3,r/R) , rf?(r) = BGi{n + 2,r/R) i?-("+2) 

\bt (!i±^) i.-(-3) , al = lBr (!i±^) , 



6 V 2 



-(n+1) 



for n > —1 . 



(37a) 
(37b) 

(37c) 

(37d) 

(37e) 
(37f) 

(37g) 



Here, C = 0.5772156649 • ■ • is Euler's constant. For a pure power-law spectrum with n < —1, 7(r) 
and (T^ diverge at long wavelength. However, Eqs. (37c) and (37d) are valid for n > —3. For 
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n < —1, c^^dj/dr = 1 in Eq. (29a). For n > — 1, cr^ converges at long wavelength but diverges at 
short wavelength with cr^ oc Rq^"'~^^^ as Rq — 0. 

The limiting behavior of ^(r) and fj{r) are given by 



fj{r) 



jrm(A)"'""' 

-nc/i(n + 2)r-("+^) , 
5i(n + 2)r-(''+^) , 



1 - 



(^^ + 0(r^) 
^^ + 0(r^) 



r<^R 
r:^ R 

r^R 
R 



(38a) 
(38b) 



For all r/R, these functions obey the relation ^ + [R^ /r)d^/dr = —nr]{r). 



The drift coefficient is plotted in Figure 5 forn G {—2, —1, 0, 1}. As expected from the previous 

section, constraints make more positive. For power-law spectra, however, o"| diverges so one 
cannot use Eqs. (31) to obtain the behavior as r ^ 0. The exact results here show that for 
—2 < n < 4, Ar ^ —oo as r ^ 0. For n < —2, ^ as r ^ oo. 

These results show that models with lots of small-scale power and substructure (n > —2 as 
k — > oo) have a strong inward drift force, while models that are smoother on small scales (n < —2) 
have negligible drift force as r — >^ 0. This is easy to understand qualitatively in the context of 
halo formation. Density fields with lots of substructure have a nonspherical mass distribution, 
with mass concentrated into subhalos. For a given halo, the mass between r/2 and r, for example, 
is concentrated into overdense substructures whose gravity field on infalling particles is stronger, 
on average, than if this mass were uniformly distributed over a spherical shell as it is for the 
average halo. Given a density correlation function that rises with decreasing r, there is a greater 
concentration of such clumps interior to r than outside of r, leading to a net inward drift force. 

Figure 6 shows the radial and tangential diffusivity for power-law spectra with n < —1. For 
n > — 1 the diffusivity diverges because the velocity dispersion diverges as A; — »• oo. For n = — 1 
the divergence is logarithmic resulting in a logarithmic dependence of Dr and Dt on radius. Figure 
6 shows that the effect of the constraints is to increase the diffusivity at radii within a few smoothing 
lengths of the peak. As we noted above, constraints on the density gradient make the diffusivities 
0. For n < — 1 the diffusivities continue to rise with increasing r because (t3^ 



negative as r 
diverges as k 

small and large k, so the asymptotic behavior is different. 



0. For the physical power spectrum used in Figures 1 and 2, converges at both 



5. Discussion of the Fokker-Planck Equation 

The use of the Fokker-Planck equation to describe dark matter clustering is novel. In order 
to gain insight into what this approach may reveal, in this section we examine the Fokker-Planck 
equation from several different perspectives. 
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5.1. Velocity Moments 

The Fokker-Planck equation describes relaxation processes in a weakly collisional fluid. In this 
regard it is similar to the Boltzmann equation, in which the right-hand side of Eq. (9) is replaced by 
a two-body collision integral that is bilinear in the one-particle distribution function. The present 
situation is physically very different however — "collisions" are not due to individTial particles but 
rather to the gravity field of density fluctuations in a Gaussian random field. Nevertheless, we 
can borrow one of the techniques used to gain insight into the Boltzmann equation, viz. velocity 
moments, to obtain a clear physical interpretation of the meaning of the Fokker-Planck equation 
and the diffusion coefficients A and D. 

Moment equations are obtained by multiplying Eq. (24) by powers of Vi and integrating over 
velocity. The three lowest moments give the mass density p, fluid velocity u, and velocity dispersion 
tensor Oij: 

p{r,t) = J (fvfif,v,t) , (39a) 
Ui (r, t) = p-^ J d^v f{r, V, t) Vi , (39b) 
{f, t) = p-^ J d^v f{r, V, t){vi - Ui)(vj - Uj) . (39c) 
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Here we are using physical (proper) spatial coordinates rather than comoving coordinates. The 
density, velocity, and velocity dispersion tensor obey the following equations: 

| + |-(P..) = 0, (40a) 

^ + «.^jn. + - — (p^,) = g^i + A, (40b) 

3 \ du ' du ' 1 d 

_ + u,—Je., + ^^e,, + ^j,, + --^ipe^,,) = 2a. . (40c) 

The first equation is the usual fluid continuity equation. The second is the usual Euler equation, 
except that the pressure is generalized to a matrix pOij because the spatial stress need not be 
diagonal for a collisionless fluid. The right-hand side has two force terms: the gravitational tidal 
field produced by p{f,t), and the drift acceleration A{f,t) caused by the correlated fluctuations 
of the mass density field. The third equation introduces a tensor 6ijk, which is defined exactly 
like 6ij except the integrand in Eq. (39c) acquires an additional factor (v^ — u^)- Eq. (40c) is a 
generalized heat equation, which requires a little more discussion. 



The velocity dispersion tensor 9ij generalizes the temperature of an isotropic gas to the case 
of an anisotropic velocity distribution. For an isotropic gas with 9ij = TSij (here T has units of v"^) 
and Dij = DSij, this equation takes the more familiar form 



3 dT w 1 ^ 
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where qi = ^pYl^=i ^ijj is the heat flux. The velocity gradient terms in Eq. (40c) generahze thepdV 
term (i.e., the T V-u term) of Eq. (41) to an anisotropic velocity distribution. For a perfect gas, the 
heat flux term is absent because the Maxwellian velocity distribution has vanishing skewness. For 
an imperfect or collisionlcss gas, however, the third moment of the velocity distribution generates a 
generalized heat flux p0ijk- Notice that the heating rate is given not by the heat flux but instead by 
its divergence. For a weakly imperfect gas, the heat flux can be approximated (e.g. in the Chapman- 
Enskog expansion approach) by a conductive flux q = —kVT where k is the conductivity, but for 
a collisionlcss gas, Oijk cannot be obtained simply from p, u, and Oij. Thus, Eqs. (40) do not form 
a closed system. This is precisely why we must use a phase space description for a collisionlcss 
gas instead of trying to modify the fluid equations. Eqs. (40) and (41), however, are pedagogically 
useful in showing us how to interpret A and D. 

The velocity diffusivity tensor Dij appears as a new term in the anisotropic heat equation, with 
no familiar counterpart in the dynamics of a collisional gas. Although Dij comes from a diffusion 
(i.e., gradient) term in phase space, it appears as a direct local heating term in real space (i.e., 
coordinate space r ). Diffusion in velocity space is equivalent to local heating: the temperature 
increases when particles spread out in velocity at a given point in space. The diffusion occurs 
in velocity space not in real space. The presence of this term begs the question: where does 
the energy come from? From the exact kinetic equation (9), it follows that the heating term is 
given by the symmetrized flrst velocity moment of the correlated force density F produced by 
fluctuating gravitational fields. From this we see that the heating comes from gravitational energy. 
For this reason, it is possible for D to be positive or negative: gravitational fluctuations can either 
locally heat or cool a gas. One should not try to draw conclusions about energy conservation 
from these results because our calculations are performed in the accelerating frame centered on 
a halo. However, the equations of motion used in deriving our kinetic theory imply local energy 
conservation. 



5.2. Langevin Equation 

The Fokker-Planck equation describes the evolution of the probability distribution for particles 
undergoing random walks. Instead of describing the system by a distribution function, we can spec- 
ify equations of motion for individual particle trajectories using the Langevin equation (Langevin 
1908). Although this is reversing the logic we used in deriving the Fokker-Planck equation, it is 
useful for providing both an understanding and a Monte Carlo method for numerically solving the 
Fokker-Planck equation. 

The Langevin equation is a stochastic differential equation for the individual particle trajec- 
tories. The Langevin equation corresponding to Eq. (24) is given by the pair of equations 

^ = v, - = 9T{r,t)+A{r,t) + Tif,t), (42) 
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where F is a force due to a zero-mean Gaussian random process with covariance 

{Ti{r, t)rjif, t')) = 2Ai(r, t) Sj^it - t') . (43) 

The velocity of a particle at time t + dt, v{t + dt), is therefore given by a deterministic term 
v{t) + [5t(^5 t) + A(r, t)]dt, plus a random term drawn from a Gaussian distribution with variance 
oc Ddt. Subtle differences exist between the cases where the variance D is evaluated at {f{t),t) 
("Ito calculus") vs. at {r{t + dt/2),t + dt) ( "Stratonovich calculus"). Details can be found in Risken 
(1989). Note that the Langevin equation requires D > 0. 

The interpretation of the Langevin equation is given by the following result (Risken 1989; 
Gardiner 2002): If a sufficiently large sample of particles is drawn at random at time to from the 
phase space distribution /(r, v, to), and their orbits are integrated using Eqs. (42), then at any later 
time t their positions and velocities are a sample from the phase space density f{f,v,t) obtained 
by integrating Eq. (24). 

Eqs. (42) are similar to the exact equations of motion for individual particles given by Eqs. (4), 
but there are several important differences. The gravity field of one realization is replaced by the 
spherically symmetric gravity field of the cnscmblc-avcragc density profile. The effects of spatially- 
varying forces in Sire represented by the drift term A and the stochastic acceleration T. These 
additional terms have the effect that, at least in the quasilinear regime (which we assumed in 
our derivation of the Fokker-Planck equation), the phase space evolution implied by the Langevin 
dynamics is equivalent to a Monte Carlo sample from the phase space density. 

This does not mean we have come full circle. For individual realizations of the ensemble, the 
density field is clustered and not spherically symmetric and the phase space is six-dimensional. 
The Langevin equations give orbits in a smooth, spherically symmetric matter distribution. The 
latter is far easier to model numerically, with much higher resolution possible, than the original 
dynamics. By taking advantage of the symmetries of the average halo we have reduced the phase 
space from six dimensions to three. 



5.3. Comparison with Classical Brownian Motion 

The Fokker-Planck equation was first written almost a century ago to describe Brownian 
motion (Fokker 1914; Planck 1917). In this case, a small macroscopic particle (a Brownian particle) 
undergoes a random walk due to its collisions with individual molecules in a liquid. Instead of the 
phase space density of molecules, the relevant quantity is the probability density W{v, t) for the 
Brownian particle to have velocity v at time t. Assuming a spatially homogeneous and isotropic 
medium, this probability distribution obeys the Fokker-Planck equation 

dW d f kBTdW\ 
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where 7 is a constant (given, for a small spherical body of mass m and radius a immersed in a fluid 

of viscosity ij, by the Stokes formula 7 = 6Trr]a/m), T is the temperature, and fee is the Boltzmann 
constant. Eq. (44) has the same form as Eq. (24) if we identify A = —jv for the drag term and 
Dij = '^{k^T /m)5ij for the diffusion term. 

We note that the Brownian motion can be equally described by the Langevin equations in 
Eqs. (42) and (43) with gx = ^, A = —^v, and Dij = ^{k-^T /ra)5ij. 

Brownian motion describes a very different physical situation than dark matter halo formation. 
Why are both described by a Fokker-Planck equation? The reason is that both problems involve 
random walks. Let us examine this more closely. In the case of Brownian motion, the random walks 
arise from impulsive collisions that cause the Brownian particle's velocity to change significantly on 
a timescale 7"^. During this time interval the particle moves a distance ^ v/j. After an elapsed 
time t, the particle undergoes 7^ steps and therefore moves a distance ~ (v'^t/j)^/'^. In thermal 
equilibrium, v"^ ~ k-QT/m so the particle moves a typical distance ^~^^/Di where D = ^k-QT/m 
is the velocity diffusivity. A random walk in space leads naturally to a diffusion equation for the 
probability distribution (Risken 1989; Gardiner 2002). 

The dark matter case at first appears to be very different. Instead of a large particle being 
buffeted by small ones, we are studying small particles being buffeted by large mass fluctuations 
(i.e. large in mass compared with the dark matter particles themselves, whose mass never enters 
our discussion). Moreover, our derivation of the Fokker-Planck equation has assumed that the 
motion is in the quasilinear regime described by Eq. (16) with a universal time-dependence for 
ip{x,t). Therefore each particle has moved in a straight line with (in appropriate units) a constant 
velocity, which is not a random walk in space. However, the velocity of a given particle is a random 
walk since tp is obtained by adding up the gravitational accelerations produced by all the mass 
fluctuations in a Gaussian random fleld. The net force on a particle is a random walk at fixed time. 

Although cosmic Gaussian random fields lead to a Fokker-Planck equation, the diffusion coef- 
ficients are very different from Brownian motion. First, the drift acceleration A is independent of 
velocity, in contrast with the viscous drag acceleration —jv for Brownian motion. In the quasilinear 
regime of gravitational clustering, the drift is not a friction at all. Instead it is radially directed 
(because of spherical symmetry for the average halo), either inward or outward depending on the 
density correlations of fluctuating substructure. Second, the velocity diffusivity tensor Dij is not 
proportional to the temperature; the dark matter is treated as being completely cold with vanishing 
temperature. In both cases the diffusivity is proportional to the mean squared displacement. In 
the cosmological case this depends simply on the power spectrum of density fluctuations and not 
on the thermal velocity of the particles. 

In the Brownian case A and D are both proportional to the Stokes drag coefficient 7. The 
proportionality between A and D is no accident but is a consequence of the fluctuation-dissipation 
theorem, which states that these coefficients have a proportionality given by the condition of thermal 
equilibrium. It is instructive to see this in a generalized case of Eq. (44) . Suppose that the Brownian 
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particle is moving in a time- independent potential $(r ). The Fokker-Planck equation then becomes 

dW _^ dW dW _ d / ksT dW\ 

dt ^ df df dv dv \ ^ m dv ) ' 

The stationary solution (for suitable boundary conditions) is given by the Boltzmann distribution 

W = iVexp[-£:/(A;Br)] , E = ^mv'^ + m$(r ) . (46) 

where A'^ is a normalization constant. With this solution not only does the left-hand side of Eq. (45) 
vanish, but so does the velocity flux on the right-hand side. Drag and diffusion balance each other 
in thermal equilibrium. Note that the mean squared velocity at fixed r equals ^k^T/m. For a 
distribution of masses of Brownian particles, thermal equilibrium results in an equipartition of 
energy with the heaviest particles moving most slowly. This is a natural consequence of having 
(D/7)^/^ equal the thermal speed. 

In thermodynamic systems, the Fokker-Planck equation governs the relaxation to thermal 
equilibrium. In the cosmological case of quasilinear clustering considered here, it is not thermal 
equilibrium but rather the statistics of Gaussian random fields that relate A and D. In Appendix B 
we investigate our cosmological Fokker-Planck equation for the relaxation to stationary solutions. 



5.4. Comparison with Globular Cluster Dynamics 

Globular cluster evolution is the most studied application of the Fokker-Planck equation in 
astrophysics. Unlike the Brownian case, gravity plays the dominant role here. A star in a globular 
cluster can be approximated as experiencing two types of gravitational forces: a smoothly varying 
potential $(r, t) due to the smoothed matter distribution in the system, and a fluctuating force 
due to many two-body interactions with other stars. The phase space density of the stars obeys 
the Fokker-Planck equation (Spitzer 1987; Binney &l Tremaine 1988) 



df ^ df df _ d 

dt df df dv dv 



This looks very similar to the cosmological kinetic equation (24) in the quasilinear regime, aside 
from a factor of two in the diffusivity that is purely a matter of differing conventions, and the 
velocity dependence of the diffusion coefficients. However, the physics of A and Dij is very different 
here. The relaxation process in globular cluster dynamics is two-body relaxation (Chandrasekhar 
1943), i.e. the dissipation arising from the fluctuating forces of discrete Newtonian gravitating 
point masses. Examining this case will shed light on the relaxation processes that can affect galaxy 
halo evolution. 

The calculation of diffusion coefficients for two-body relaxation is described in § 8.3 of Binney 
&; Tremaine (1988). To understand the results it is helpful to consider the process of Chandrasekhar 
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dynamical friction acting on a test particle with mass rrit and velocity vt moving through a spatially 
homogeneous sea of background particles of mass rrif, and velocity distribution f{vb) with mass 
density pb = nib f f{vb)d?Vb. For an isotropic Maxwellian distribution f{vb) oc exp(— i;^/2(T^), the 
test body feels an acceleration 



dvt 
dt 



A-KG'^Pbirnt + rrib) In A 



9 X 

erf(X) - =^e- 

TT 



Vt 



(48) 



where X = vt/ {\/2ab)- 



Now treat each star within the globular cluster one at a time as a test particle, with all the 
other stars serving as background particles. If the velocity distribution is Maxwellian and spatial 
gradients are ignored, dynamical friction will contribute A = —dv/dt given by Eq. (48) with vt 
replaced by v. More generally, for a spatially homogeneous distribution on scales larger than the 
interparticle spacing, two-body relaxation gives (Binney &: Tremaine 1988) 



d 

A = AnG'^mbimt + m^) In A -wzHv ) 



D. 



4.GV,lnA ^^^^^ 



dv 
■9{v), 



(49a) 
(49b) 



where In A is the Coulomb logarithm, and the Rosenbluth potentials h and g are given by (Rosen- 
bluth, MacDonald, & Judd 1957) 



hiv) ^ 



f{vb)d^Vb 

\v -Vb\ 



9{v) = J fivb)\v-Vb\d^Vb. 



(50a) 
(50b) 



For the isotropic Maxwellian distribution f{vb) oc 6xp(— t;^/2c7^), A = —j{v)v is a drag force and 
the needed coefficients are 
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where X is defined as before (with vt ^ v), and Dij is decomposed into two components: 



(51a) 
(51b) 
(51c) 



A.- = A^Dii + - {5ij - -^J . (52) 



We note that the drift A contains two terms: the first term is proportional to the test particle 
mass mt while the second is independent of mj. The first term represents a "polarization cloud" 
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effect, which arises from the fact that a test particle moving relative to a background deflects more 
upstream than downstream background particles (Gilbert 1968; Mulder 1983). This process results 
in a density gradient and hence a drag force on the test particle, whose amplitude increases linearly 
with nit since the gravitational deflection is proportional to nit- In contrast, the second term in 
A and the diffusivity tensor D in Eq. (51) are both oc Pbmi, and independent of mj. These two 
terms are due to fluctuations in the gravitational potential caused by granularity in the background 
particle distribution. 

The derivation of the Fokker-Planck diffusion coefficients for two-body relaxation is completely 
different from our derivation based on quasilinear cosmological density fluctuations. In both cases, 
fluctuating forces lead to random walks in velocity space, hence to the Fokker-Planck equation. 
The nature of the fluctuating forces and their consequences, however, are very different in the two 
cases. 

Our cosmological Fokker-Planck derivation is exact (in second-order cosmological perturbation 
theory); it yields no drag; and the resulting drift and diffusion coefficients depend on position but 
are independent of velocity, as shown in Eqs. (28)-(30). In contrast, the two-body relaxation 
calculation is based a Taylor series expansion of the phenomenological (not exact) master equation; 
it yields no radial drift; and the drag and diffusion coeflficients depend on velocity but are (in the 
usual derivation) independent of position (Rosenbluth, MacDonald, &; Judd 1957). Furthermore, 
the two-body relaxation rates are proportional to while the rates we compute are proportional 
to G from Eq. (28). In the cosmology case, the phase space correlations are built into the initial 
power spectrum of density fluctuations, whose amplitude is given independently of G. With two- 
body relaxation, however, correlations are gravitationally induced by fluctuations instead of being 
present ab initio, hence the rates pick up another factor of G. 

The differences in the A term (radial drift versus drag) have an important consequence for 
equipartition of energy. As shown previously in the discussion of Brownian motion, drag and 
diffusion work together to drive a system to thermal equilibrium in which the mean squared velocity 
is proportional to kT/m. Binney and Tremaine show (p. 513 of Binney & Trcmaine 1988) show that 
this also happens with two-body relaxation. Because equilibrium systems with only gravitational 
forces generally have negative specific heat, there is no stable thermal equilibrium. Nonetheless, 
two-body relaxation tends to drive the velocity distribution toward the Maxwellian form, with 
equipartition of different mass species. 

In the quasilinear cosmological case, the drift and diffusion coefficients are independent of both 
velocity and particle mass. Therefore the equilibrium velocity distribution cannot depend on parti- 
cle mass — there is no equipartition of energy. This result is similar to violent relaxation, a process 
arising during the initial collapse and virialization of a dark matter halo when the gravitational po- 
tential is rapidly varying in time (Lynden-Bell 1967). Our relaxation process, however, arises in the 
quasilinear regime and is present even if d(f)/dt = as it is (at fixed comoving position) if fim = 1. 
Moreover, violent relaxation (and its partner, phase mixing) are generally understood as arising 
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from the collisionless dynamics represented by the left-hand side of Eq. (47) , while both two-body 
relaxation and our cosmological relaxation process are represented by the "collision" terms on the 
right-hand side. 

Why did our cosmological calculation not yield a drag term (and a corresponding diffusion 
term linked to it by the fluctuation-dissipation theorem)? We suspect that it is because our calcu- 
lation was limited to small-amplitude perturbations about a homogeneous and isotropic expanding 
cosmological model. The induced effects of substructure would only come in at higher order in 
perturbation theory. In the nonlinear regime we expect two types of drift terms to be present: 
A = af — ^v, where a is the radial drift and 7 is the drag coefficient. The Chandrasekhar calcula- 
tion suggests that the drag and its accompanying diffusivity will depend on both position (through 
Pb) and velocity. The radial drift is absent in the Chandrasekhar calculation because the back- 
ground there is assumed to be spatially homogeneous; with an inhomogeneous distribution in a 
virialized halo there should be a radial drift (and corresponding diffusivity) that could depend on 
both position and velocity. In the nonlinear regime drag should arise from dynamical friction just 
as with the globular clusters, except that the "background" particles whose discreteness causes the 
TTife terms in Eqs. (51) are now the subhalos and substructures that rain upon a halo. Because 
we have not yet extended our derivation to the fully nonlinear regime, these expectations, while 
plausible, have yet to be demonstrated. 



6. Summary and Conclusions 

In this paper we have developed a cosmological kinetic theory, valid to second-order in pertur- 
bation theory, to describe the evolution of the phase-space distribution of dark matter particles in 
galaxy halos. This theory introduces a new way to model the early phases of galaxy halo formation, 
which has traditionally been studied by analytic infall models or numerical N-body methods. 

The key physical ingredients behind our kinetic description are stochastic fluctuations and 
dissipation caused by substructures arising from a spectrum of cosmological density perturbations. 
The kinetic equation that we have obtained at the end of our derivation, Eqs. (24)-(26), has 
the standard form of a Fokker-Planck equation. The diffusion coefficients A and D represent 
acceleration due to a drift force and velocity-space diffusion, respectively. To second order in 
perturbation theory, these coefficients are related to various covariance matrices of the cosmological 
density and velocity fields given by Eqs. (25)-(26). These expressions can in turn be written in 
terms of the familiar linear power spectrum P{k) of matter fluctuations, as shown in Eqs. (29) and 
(30). 

The results for A and D from second-order perturbation theory are shown in Figures 1-4 for the 
currently favored ACDM model, and in Figures 5 and 6 for various scale-free models with power-law 
P{k) oc A;". Our results indicate that dissipative processes are important during the initial collapse 
of a dark matter halo, and it is not valid to model halos using the spherically symmetric Vlasov 
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(coUisionless Boltzmann) equation. Furthermore, we find that the diffusivity is initially negative 
close to the center of a halo, indicating a thermodynamic instability. The source of this instability 
is gravity: perturbations enhance gravitational collapse of dark matter halos. 

Wc emphasize that our derivation leading to the Fokkcr-Planck equation given by Eqs. (24)- 
(26) is exact to second order in cosmological perturbation theory. Wc do not follow the frequently 
used approach (e.g. Binney &: Tremaine 1988) that begins by expanding the collisional terms 
of the non-exact master equation in an infinite series of changes in phase-space coordinates (i.e. 
A.w), and arrives at the Fokker-Planck equation by assuming weak encounters (i.e. small \Aw\) 
and truncating the series after the second-order terms. Instead, our starting point is Eq. (9), the 
first BBGKY hierarchy equation for a smooth one-particle distribution function. This is an exact 
kinetic equation; solving it requires specifying the one-particle distribution at an initial time and 
the two-particle correlation function at all times. We have shown in § 3 and Appendix A that if the 
cosmological density field 6 and displacement (or velocity) field i/j are Gaussian random fields, the 
resulting two-particle correlation function depends on the one-particle distribution in a simple way 
given by Eq. (A22). Eq. (A22) is exact to second order in perturbation theory and provides the 
closure relation for the BBGKY hierarchy in Eq. (9). Combining Eq. (A22) with Eq. (9) results in 
Eq. (24), which is a Fokker-Planck equation. 

The cosmological dissipative processes identified in this paper are quite different from the 
standard two-body relaxation, which also leads to a Fokker-Planck equation. We find that there is 
no dynamical friction (i.e. terms proportional to —v in the diffusion coefficients) in second-order 
cosmological perturbation theory. Instead, there is a radial drift toward (or away from) the center of 
a halo due to the clustering of matter within the halo. The usual treatment of two-body relaxation 
has no such drift term because it assumes the medium to be homogeneous. Thus we have identified 
a new relaxation process that must affect the phase-space structure of dark matter halos. 

Although our derivation is exact to second order in the density fiuctuations, it is valid only 

when the flTictTiations are small. We are therefore describing only the early stages of dark matter 
halo formation. The expressions we obtain for the diffusion coefficients are plainly invalid in the 
nonlinear regime, where we expect dynamical friction to be present as well as radial drift. In a later 
paper we will investigate the nonlinear generalizations of the exact second-order results derived in 
this paper. 

The Fokker-Planck description should still apply in the nonlinear regime provided that the 
fluctuating force on a particle can be modelled as a Markov process consisting of a random walk 
of many steps. We conjecture that this description will be approximately valid when the matter 
distribution is modelled as a set of clumps (i.e., the halo model) which scatter individual dark 
matter particles away from the orbits they would have in a smooth, spherical potential. The 
Fokker-Planck description should apply not only to the initial collapse and virialization of a dark 
matter halo containing substructure, but also to the subsequent evolution under minor mergers. 
Our inability to solve exactly the nonlinear dynamics, however, means that the diffusion coefficients 
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will have to be calibrated using the nonhnear halo model or numerical N-body simulations. 

Calibrating the diffusion coefficients is not equivalent to reproducing the N-body simulations. 
N-body simulations provide a Monte Carlo solution of the BBGKY hierarchy, which is much more 
general than the Fokker-Planck equation. If the simulation results are describablc by Fokkcr-PIanck 
evolution, this already represents a significant new result. Moreover, even if functional forms for the 
radial dependence of the diffusion coefficients have to be calibrated with simulations, the solution of 
the Fokker-Planck equation contains far more information because it gives the complete phase-space 
density distribution, not simply radial profiles. 

Another important step will therefore be to actually solve the Fokker-Planck equation. In 
deriving it we have identified and quantified the relaxation processes affecting halo formation and 
evolution, but this is only the first step. In Appendix B we presented solutions in unrealistically 
simple cases just to highlight the roles played by radial drift, drag, and diffusion. An outstanding 
question is whether these processes actually arc rapid enough to erase the memory of cosmological 
initial conditions sufficiently so that dark matter halos relax to a universal profile. Once we have 
a nonlinear model for the diffusion coefficients, we will address this important question. 

The work presented in this paper has a number of implications. First, wc have highlighted the 
importance of the phase space density for understanding dark matter dynamics, and have developed 
a method for calculating its evolution during the early phases of structure formation. Second, 
we have identified a dissipative process that may erase memories of initial conditions during the 
early phases of galaxy halo evolution. This erasure is necessary for universal halo density profiles, 
although further investigation is needed to quantify how this dissipative process affects the density 
profiles. Third, the statistical description of the potential fluctuations caused by substructure 
that we have developed can also be added to investigate related problems such as the heating of 
galactic disks (Benson et al. 2003) and the merging of central black holes when galaxy halos merge 
(Kauffmann & Haehnelt 2000; Haehnelt & Kauffmann 2002; Hughes & Blandford 2003). 
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from the Research Corporation, and NASA grant NAG5-12173. 

A. Statistics of Constrained Gaussian Random Fields 

In this Appendix we derive expressions for the probability distributions of displacement and 
density perturbation, p{ijji) and p{d2,il^i), that are needed in Eqs. (21) and (22). We then use 
them to evaluate Eqs. (21) and (22) at high redshift when the matter distribution was described 
by a Gaussian random field of small-amplitude density perturbations. We use comoving spatial 
coordinates x = f/a{t). 
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In accordance with the standard cosmological model, we assume that S{x,t) is a Gaussian 
random field with power spectrum P{k,t), and that the displacement and density are related by 
Eq. (17). These arc the appropriate assumptions at high redshift in standard models of structure 
formation. Strictly speaking, at high redshift, when linear theory applies, the evolution of the dis- 
tribution function is trivial because the density and velocity fields evolve by spatially-homogeneous 
amplification. Here we use the statistics of Gaussian random fields to compute quantities driving 
the lowest-order corrections to linear evolution. We will find that the fluctuating force is formally 
second order in perturbation theory. 

A Gaussian random field is described most simply by its Fourier transform, which we define 

by 




We suppress the time dependence for convenience. Before we apply a constraint, both S{x) and 
S{k) have zero mean. The covariance of d{k) gives the power spectrum, 

{S{ki)5ik2)) = (27r)3p(fei)5D(fei + k) , (A2) 

where Sd is the Dirac delta function. In the linear regime, ^{x ) is a linear functional of d{x). 
The multivariate distribution of any linear functional of a Gaussian random field is itself Gaussian 
(i.e. multivariate normal) and is therefore described completely by its mean and covariance matrix. 
This remains true even for a constrained Gaussian random field (Bardeen et al. 1986). 

In our case, we wish to constrain the initial density field to represent a proto-halo. The 
constraints are applied to the smoothed density field 

A(f ) = J Wr{x - x' )5{x' ) = j ^ WR{k)6{ky'^-^ , (A3) 

where 

j^WR{k)e''-^ (A4) 

is a smoothing window with unit volume integral, f d'^x Wr{x ) = 1. The subscript R indicates the 
smoothing length. We assume that the smoothing window is spherically symmetric so that WR(k) 
is real and depends only on the magnitude of the wavevector. Throughout the paper we will use 
the Gaussian window Wr{x) = eK^{-x'^ /2B?) and WR{k) = ex.Y>{-k'^B? /2). 

Unfortunately there is no analytic theory to tell us where halos will form. However, there are 
a variety of plausible approaches. For example, following Bardeen et al. (1986) we might constrain 
A to have a maximum of specified height, orientation, and shape at x = 0. This peak constraint 

involves ten variables (A and the components of its first and second derivatives with respect to 
each of the coordinates), presenting a formidable challenge to evaluating the covariance matrix of 
the four variables {52-,'4^i)- The peak constraint is appealing but it is difficult to calculate and in 
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practice it has not been found to predict well the actual formation sites of dark matter halos (Katz, 
Quinn, & Gelb 1993). 

Instead of trying to impose complicated constraints for halo formation (e.g. Monaco et al. 
2002), we err on the side of simplicity by choosing to use only the zeroth and first derivatives of 
the smoothed initial density field A(x). We require that 5; = be an extremum of the smoothed 
density field by requiring VA(0) = 0. Of course, an extremum may be a maximum, a minimum, 
or a saddle point. We therefore keep track also of the smoothed density Aq = A(0). Provided that 
Aq is greater than 2 standard deviations, the extrema of A(i;) are close to the peaks (Bardeen et 
al. 1986). 

How will we know if our results depend strongly on the choice of constraint? We will do this 
in two ways. First, we will vary the height of the extremum, Aq. Second, we will drop the gradient 

constraint and use the single constraint Aq. While these tests will not prove that we have a good 
model for the sites of halo formation, they do provide control cases for us to assess the sensitivity 
of our results to the details of the initial constraints. 

Now we must calculate the mean and covariances of (^25'0i) subject to constraints (Ao,Ao) 
where 

Ao = A(x = 0), Ao;eVA(x = 0). (A5) 

The method for this calculation is based on the theorem presented in Appendix D of Bardeen et al. 
(1986), which states that if Yj^ and Yb are zero-mean Gaussian variables (more generally, vectors 
of any length), then the conditional distribution for Yb given Ya, p{Yb\Ya) = p{Ya,Yb)/p{Ya), is 
Gaussian with mean 

{Yb\Ya) = J dYB p{Yb\Ya) Yb = {Yb ® Ya) {Ya ® ^a)"' Ya (A6) 
and covariance matrix 

C{Yb,Yb) = j dYBp{YB\YA) AFs ® AF^ = {Yb ® Yb) - {Yb ® Ya) {Ya ® Ya)"' {Ya ® Yb) , 

(A7) 

where AYg = Yb — {Yb\Ya)- The symbol ® denotes a tensor product and angle brackets denote 
mean values. Thus, (Ya^Ya) is the covariance matrix of the constraints while C(Yb,Yb) = 
(AYb (g) AYb)c = (AYb (8) AYb| Ya) is the covariance matrix of Yb subject to the constraints. (The 
symbol C, used either as a function or as a subscript, implies that a constraint is applied.) The 
juxtaposition of two matrices implies matrix multiplication. 
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A.l. Covariance of Constraints and Variables 

To apply these results, first we compute the covariance matrix of constraints Ya = (Aq, Aq). 
Using Eqs. (A1)-(A4), we get the covariances 

^P{k)WUk), (A8a) 
(AoAo) = 0, (A8b) 
(Ao Ao) = all, aj^-j k^P{k)Wl{k) . (A8c) 

where I is the unit tensor. 

Next we consider the covariance matrix of the variables Yb = ^2) V^i)- Since 5 = —{d/dx)-^, 
the displacement field ip has an additional factor ik/k"^ in the integrand of Eq. (Al). For the off- 
diagonal elements of the covariance matrix of Yb we obtain 

^ P(fc)io {kr) , r = Xi-X2, ( A9a) 

5iV7i) = 0, (A9b) 
^i) = -V{r)r, V{r) ^ I ^ Pik)'^ , (A9c) 

where the spherical Bessel functions are jo{x) = a;~^sinx and ji{x) = a;~^(sina; — xcosx) = 
—djo/dx. The diagonal variances are 

{Si) = {si) = m , (AlOa) 
^ ^ \ 1 f d^k 

V'l ^ V'l) = 41 ' 4 = 37 k-^P{k) . (AlOb) 

Note that in Eqs. (A9) (and in the following) we are using r as a comoving displacement vector and 
not as the radius vector in proper coordinates. Proper coordinates are not used in this Appendix. 
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A.2. Cross-covariance between Variables and Constraints 

Finally, we have the cross-covariances between our primary variables and the constraints, i.e. 
the components of (Yb <8) Ya)- These follow from Eqs. (A1)-(A3): 

^P{k)WR{k)jo{kr) , r = \x\, (Alia) 

^{x)Ao) = -n{r)S, fj{r) ^ I -0^ P{k)WR{k)^-^ , (Allb) 

S(x)Ao) = -^r , f = x/r, (Allc) 
/ dr 

if(x)®Ao) = r^f^r + fjl. (Alld) 
/ dr 

To arrive at the last expression above, we have used the identity 

'jiikr) 



I 



— e ^^^n®n = r— 
47r or 



kr 



r®r + 4^I, (A12) 
kr 



where k = kn, x = rr and I is the unit tensor. 



A. 3. Mean of Variables Subject to Constraints 

Using Eqs. (A6) and (A8)-(A11), we get the mean values subject to the extremum constraints 
(Ao,Ao = 0), 

{5i)c = ((5i|Ao,Ao) = %"(ri) , n = , (Al3a) 
{S2)c = (52|Ao,Ao) = ^^>2) , r2 = \x2\, (Al3b) 
{ifi)c = Ui\Ao,Ao) = -^fj{ri)xi . (A13c) 

The subscript C denotes "subject to constraints." 

Eqs. (A13) would be unchanged had we dropped the extremum condition and imposed instead 
the single constraint A(0) = Aq; the constraint Aq = has no effect on the mean values. This 
follows because (AqAq) = so that the only way that Aq can enter the mean values is through 
terms proportional to Ao/af which vanish when the extremum constraint is imposed. 



A. 4. Covariance of Variables Subject to Constraints 

We now compute the covariance matrix of our variables Yb = ((^i, (^2, V'l) subject to the ex- 
tremum constraints = (Aq, Aq). The covariances follow from Eqs. (A7)-(A11). Recalling the 
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definition C{Yb,Yb) = (AYb AYb\Ya), we obtain 



C{5i,d2) = e.(r) 2 2 3— 3— ''I ■ ^2 , r = xi- X2 , ri = — , r2 = — , (A14a) 

erf ari ar2 ri r2 



C((52,V^i) = -r/(r)f+^^ + ^-^ [(C"i-3f?i)(fi-f2)fi + f?if2] , (Al4b) 

(Tq (tJ ar2 

C((5i,V^i) = fl^ + i^^(ei-2f7i)xi , (A14c) 
C(V'i,V'i) = c^ri (g) fi + ct(I - fi (g) fi) , 

nr/A^ /Ci - 2r7i 



. .^-(^^j -^^j . (A14d) 
ct(n) = "l-i^f^ • (Alio) 

where = C('^i) aiid ryi = 77(ri), and we have used the relation ri{df]i / dri) = ^1 — St^i- Note 
that C((5i,V'i) is a vector and C(V'i,^i) is a second-rank tensor. Later we will need the tensor 
C~-^('0i, ^1), which is the matrix inverse of C{i/ji,iI^i)- 

The role of constraints is easy to see in Eqs. (A14). The constraint on Aq is responsible for 
the terms proportional to a^"^ while the constraint on Aq gives rise to the terms proportional to 
crf^. Therefore, the constraints are easily removed by setting crj"^ — > (to eliminate the constraint 
on Aq) or a^^ (to eliminate the constraint on Aq). As expected, in the absence of any constraints, 
the covariance matrix recovers the same expressions as the unconstrained covariance in Eqs. (A9) 
and (AlO). 

We will need one more covariance, C{ipi,ipo) subject to the extremum constraints where tpQ = 
if{0). Using Eqs. (A7), (All), and (A12), we get 

<^(V'i,V'o) = T^n®ri + — (I-n(8)ri) , j{ri) = P{k)k-^ ji{kri) . (A15) 

dri ri J {2Try 



A. 5. Probability Distributions 

Now we are ready to calculate the probability distribution functions and means of various 
quantities needed in § 3. Specifically, we needp(V'i), {Si\'4'i)cj {^2\'4'i)cj {^2)0^ and ((Ji52|V'i)c that 
appear in the key Eqs. (21) and (22) for f{wi,t) and f2c{'Wi,X2,t). Eq. (A13b) already gives {62)0 
but we need to calculate the other quantities. 

First, the probability density distribution for the displacement p{ipi) (subject to the extremum 
constraints) is a three-dimensional Gaussian whose mean and covariance are given above, implying 



-1/2 



exp 



(AI6) 
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Note that the argument of the exponential is a scalar — the dot denotes a contraction (dot product) . 

The joint distribution of displacement and density (subject to the constraints) is a little more 
complicated to work out since this requires inverting the 4x4 covariance matrix of {S2,ipi)- The 
distribution factors into conditional and marginal distributions, p{S2,ipi) = p{ipi)p{S2\ipi) ■ Straight- 
forward algebra gives the conditional distribution, 



p(<52|V^i) = (27rQ)-i/2gxp 



1 

2Q 



(A17) 
where 

Q = C{S2,62) - C{S2, ^i) ■ C-^ifu^fi) ■ C{tPu 62) . (A18) 

Note that the second term in Q (and the similar term in eq. A17) is a scalar. Writing it out with 
indices and using the summation convention, €{62, 62) — Q = Ci{S2,tpi)C^^{tl)i,il'i)Cj{ilJi,62)- 

Eqs. (A17) and (A18) could also have been derived by noting that {Yb\Ya)c = {^2\'^i) c — {^2} c 
for Ya = '01 — {tpi)c and ¥3 = 62 — {S2)c- Then using Eq. (A6) and relations {YbYa)c = 
{{S2 - {52)c){i^i - ('0i)c))c = C((52,V'i) and {YaYa)c = C('0i,V'i), we obtain 



{h\'$i)c = {S2)c + C{62,tPi)-C-\if,Ji)-{A-{A) 



c) 



= (52)c-C(52,V^i)--^logp(V^i) , (A19) 

dtpi 

which is identical to the mean in Eq. (A17). The constrained mean (5i|V'i)c follows simply by 
replacing {62)0 with {Si)c and (7(^2, V'l) with C{Si,ipi). 

The only remaining constrained mean that we need to evaluate is {Si62\'4'i) c ■ Direct evaluation 
of the conditional distribution would require inverting a 5 x 5 matrix. It is much simpler to use 
Eqs. (A6) and (A7) with Yb = {{Si - {Si)c), {S2 - {S2)c)} and Ya = tfi - {ipi)c- Prom the 
off-diagonal element of (AYg (g) AYb|Ya), we get 

{hh\i^i)c = {Si)c{h)c + C{6i,52) - C{5iM) ■ C-\^fiM) ■ Ci^fl,62) . (A20) 



At last we are ready to evaluate the one- and two-particle distribution functions for small- 
amplitude Gaussian fluctuations. Substituting ((5i|'0i)c into Eq. 21), to second order in S we 
obtain 

f{w, t) = pil + {d)c) - C{S, ^ ) ) (Hb)-^ + 0{6') , (A21) 

where {S)c is given by Eq. (Al3a) with ri = r and C {6, ip) is a vector equal to the covariance of the 
constrained 6 and ■0 given in Eq. (A14c). Without an initial constraint we would have {5)c = and 
C{S,'ip) = 0. An initial constraint on the proto-halo shifts the mean density and causes the density 
and velocity to be correlated hence it changes the phase space density. Note that Eq. (A21) implies 
that in linear theory the spatial and velocity distributions decouple. As expected, the mean spatial 
density is p{l + {S)c)- The velocity distribution at each point in the halo is simply the underlying 
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Gaussian of Eq. (A16) with the mean shifted by C{S,ip). Note that the value of the proto-halo 

smoothed density constraint Aq affects the density profile, {S)c oc Aq, but it does not affect the 
velocity distribution up to second order in perturbation theory. The imposition of the extremum 
constraint Aq = affects the phase space density only through the term in C{6,tp) (x aj'^. Setting 
crj"^ —>■ everywhere removes the constraint on Aq. 

The two-particle correlation is found by substituting Eqs. (A19) and (A20) into Eq. (22). The 
result is 

f2c{wi,X2,t) = p^{c{di,S2)-C{6i,^i)-C-\^i,^i)-C{^i,52) 

- [c(<52,V^i)-C(5i,V^i)(<52)cl ~logp{^fl)]p{i^l)iHb)-^ . (A22) 
This equation is the main result of this Appendix. It is exact for Gaussian random fields. 

B. Stationary Solutions and Relaxation to Equilibrium 

Stationary solutions of the Fokker-Planck equation are relevant for systems that relax on a 
timescale faster than external changes take place. The stationary solution for classical Brownian 
motion was given in § 5.3. Stationary solutions may also exist for cosmological systems where the 
constant-temperature heat bath is replaced by small-scale structure that rains upon dark matter 
halos. 

Stationary solutions have df/dt = and F = Af — Y) ■ df /dv = so that the fiux vanishes 
in velocity space. We consider spherically symmetric solutions as appropriate for the average halo. 
Prom Jeans' Theorem, the distribution function is a function of the integrals of motion (Binney & 
Tremaine 1988). For spherical systems it is often assumed that stationary solutions have the form 
f = f{E,L) where 

E = ^v'^ + ^{f) , L=\rxv\. (Bl) 

(The particle masses are irrelevant, so we use the specific energy and angular momentum.) The 
collisionless Boltzmann equation gives no constraint on f{E,L). The Fokker-Planck equation, 
however, constrains / because stationarity requires that D • dlnf/dv = A. Assuming spherical 
symmetry, we write 

A = af - , (B2) 

introducing the radial drift coefficient a and drag coefficient 7. We assume isotropic diffusivity 
D. For the globular cluster case with Chandrasekhar's dynamical friction, a = 0. For the case of 
quasilinear cosmological fluctuations arising from Gaussian random flelds, we have seen that 7 = 0. 
(The cosmological case also allows for different diffusivities in the radial and tangential directions, 
but that is an unnecessary complication here.) 
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Requiring the phase space flux to vanish now gives the fohowing constraints on the distribution 
function f{E, L): 

dlnf 7 a r^dlnf a 

It is interesting to see what conditions lead to thermal equilibrium with a Maxwellian distribu- 
tion of velocities, i.e. Eq. (46) with / = W. Eq. (B3) shows that the conditions are a = 

and j/D = m/{k^T). Thus, the drift must be velocity dependent and behave like a drag, 
A cc —V. Friction converts shear into heat in such a way as to drive the system towards ther- 
mal equilibrium. Chandrasekhar dynamical friction satisfies these conditions; therefore, two- 
body relaxation drives the velocity distribution towards the Maxwell-Boltzmann form. Quasi- 
linear cosmological fluctuations, however, are quite different, having 7 = and a 7^ 0, implying 
dlnf/dE = —{r^/L)dlnf/dL = a/{vrD). While these equations do not have simple general solu- 
tions, it is clear that cosmological fluctuations do not drive a system towards thermal equilibrium 
with equipartition of energies. Our cosmological relaxation process is similar to violent relaxation 
in this respect (Lynden-Bell 1967). 

Although the general cosmological case of relaxation in a dark matter halo with a 7^ is 
complicated and cannot be fully solved here, there is a special case worthy of closer examination, 
namely, the evolution of an unconstrained cosmological Gaussian random field, where fluctuations 
generate peculiar velocities. That is, we drop the constraint of having a halo centered at r = 
and instead consider the velocity distribution of dark matter particles at a randomly selected point 
in space. To consider the cosmological case, wc should use the appropriate comoving variables to 
factor out the Hubble expansion. In place of f and v, we use the following variables: 

X = a~^f , u = a{v — Hr) , (B4) 

where a(t) is the cosmic scale factor and H = dlna/dt. With this change of variables, the Fokker- 
Planck equation (24) becomes 

If the velocities are measured in the comoving frame, is replaced by g. When using comoving 
coordinates the gravitational field is reduced by the Hubble acceleration, {d^a/dt^)x. If the mean 
density field is homogeneous, then all particles move with the unperturbed Hubble expansion, 
r oc a{t), implying g = {d^a/dt^)x. Thus, both the velocity and acceleration terms on the left-hand 
side vanish when df/dx = 0. If we apply no initial constraint on the density field, then the one-point 
distribution function / is a function only of time and peculiar velocity v — Hr. We use the variable 
u for the velocity because v — Hr decreases in proportion to (t) for a freely-falling body in 
an unperturbed homogeneous and isotropic expanding universe. If not for the right-hand side, the 
solution of Eq. (B5) for a homogeneous and isotropic model would be f{x, u, t) = f{x = 0, u, t = 0), 
i.e. the the initial peculiar velocity distribution would simply redshift with cosmic expansion and 
remain spatially homogeneous. 
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The Fokker-Planck equation does not, however, describe an unperturbed homogeneous and 
isotropic universe. Rather it describes the statistical average of an ensemble of universes with 
fluctuations. In the case of unconstrained initial conditions, / = f{u,t) and Eq. (B5) reduces to 



where = and D = D{u,t) in general. We have used statistical isotropy to require 

A = —a^u where u is a unit vector in the direction of u. In second order perturbation theory 
(§ 4.1), au = (no drag) and D = ATiGpaHba'^ (assuming that we measure velocities in the 
comoving frame.) The Fokker-Planck equation then takes the form of a diffusion equation with D 
independent of u, and the general solution is 



The effect of fluctuations is simply to spread the initial velocity distribution (if D > 0) while 
retaining the Gaussian form. Physically, the gravitational fluctuations of substructure cause local 
heating. 

Contrary to Eq. (B7), cosmological simulations show that the distribution of peculiar velocities 
in the nonlinear regime is closer to exponential rather than Gaussian, at least in the tails (Sheth 
Sz Diaferio 2001). It is interesting to ask whether the Fokker-Planck equation can explain such a 
result. With the drift term retained, Eq. (B6) has the stationary solution 



If au/D is positive and independent of n, the result is indeed an exponential distribution of peculiar 
velocities u. In second-order perturbation theory we found that = and D is independent of u, 

however Chandrasekhar's dynamical friction gives A/D oc u. Evidently, a third case is needed to 
explain the exponential velocity distribution — there must be a drag such that Uu/D is independent 
of u, quite unlike dynamical friction from two-body relaxation. Determining whether such a drag 
term arises from nonlinear gravitational clustering is a subject for further work. 

As one final example, we consider the effects of both drag and radial drift in a simplified 
model illustrating the effects of collisions on a system for which the relaxation time is much longer 
than the dynamical time (e.g. the crossing-time for a particle in a halo). Under these conditions, 
the system at all times is nearly in coUisionless equilibrium with the phase space density being a 
function of the actions alone (and not the angles, using action-angle variables). Because the actions 
are conserved in the absence of collisions, for quasi-static evolution Eq. (24) simplifies to 




(B6) 




(B7) 




(B8) 
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where we have generalized the radial drift to an arbitrary vector field a which we assume to be 
independent of v. We can solve this equation exactly for the simple case where a, 7, and D are all 
constants. The general solution is the Green's function solution 



f{r, v,t) = J <fv' G{r, v - v' , t) f{r, v',0), 



with Green's function 



G{f,v,t) 



(27r 



-g exp < IS 



(1 



a 



27 



s-D-s(l-e-2^*)| 



(BIO) 



(Bll) 



Gif,v,t) 



(B12) 



The parameter 7 ^ is a relaxation time. Eq. (Bll) simplifies to a normalized Gaussian in two 
limits, 

jA/'(at,2Dt) , 7t<l, 
\Ar(a/7,D/7), 7i»l, 

where M{u, M) is a multivariate normal with mean u and covariance matrix M. The solution for 

(3 = 7 = was given previously in Eq. (B7). Now wc have a more general solution showing that 
the equilibrium (for constant diffusion coefficients) is a MaxwcUian with velocity dispersion tensor 
D/7 just in Eq. (B8) if a„ = 'ju, but now the drift vector a produces an offset of the mean velocity 
(v) = a/'j. Thus, velocity-independent drift a and velocity-dependent drag —'jv play completely 
different roles in collisional relaxation. In later papers we will investigate these roles for more 
realistic models of nonlinear halo evolution. 
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Fig. 1. — Radial drift coefficient plotted as a function of halo radius r for Gaus- 
sian smoothing lengths i? = 1, 0.5, 0.25, 0.1, and 0.05 Mpc. It is computed from 
Eq. (29a) for the ACDM model with (O^, J^a, O;,, /i) = (0.3,0.7,0.05,0.7). For compar- 
ison, the dashed curve shows Ar for an unconstrained field. The vertical axis shows the 
dimensionless A^/a^, i.e. the drift normalized by the rms of the peculiar gravitational 
acceleration. 
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Fig. 2. — Diffusion coefficient D in the radial (r) and tangential (t) directions plotted 
as a function of lialo radius r for Gaussian smootliing Icngtlis i? = 1, 0.5, and 0.05 
Mpc. It is computed from Eqs. (29b) and (29c) for the same cosmological model as 
Fig. 1. The dashed curve in each panel is for an unconstrained field. The vertical axis 
shows the dimensionless Dij/a'^, i.e. the difFusivity normalized by the mean squared 
peculiar gravitational acceleration. 
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Fig. 3. — The same as Fig. 1 except that no extremum constraint (i.e., Aq = VA(x = 
0) = 0) is apphed. 
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Fig. 4.— 



The same as Fig. 2 except that no extremum constraint is apphed. 
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Fig. 5. — Radial drift coefficient Ar{r) for power-law spectra, P{k) oc A;", plotted as a 
function of halo radius r scaled to the Gaussian smoothing radius R. Results are shown 
for four spectral indices, n = —2, —1,0, and 1, and for constrained and unconstrained 
fields. The amplitude of Ar corresponds to P{k) = 27r^i?/c" normalized to B = 1. 
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Fig. 6. — Diffusion coefficient D in tlie radial (r) and tangential (t) directions for 
power-law spectra, P{k) oc A;", plotted as a function of halo radius r scaled to the 
Gaussian smoothing radius R. Results are shown for two spectral indices, n = — 1 and 
—2, and for constrained and unconstrained fields. For n < — 1, constraints cause the 
diffusivity to approach a negative constant for r <^ R. The amplitude of D corresponds 
to P{k) = 27r^5/c" normalized to B = 1. 
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